Does the air pollution of a country have a direct influence on the annual average temperature?¶
1. Business Understanding:¶
We all know that air pollution is a factor in climate change. What we also know is that this is a global problem. But could it also be that the emissions of a country could already have a direct influence on the annual average temperature of that country? With data from the Community Emissions Data System (CEDS) and temperature data from FAOSTAT Temperature Change, this statement will be verified within a CRISP-DM analysis.
CEDS: https://github.com/JGCRI/CEDS
FAOSTAT: https://www.fao.org/faostat
Data Sources
CEDS-Data: https://ourworldindata.org/explorers/air-pollution
FAOSTAT-Temperature-Data: https://www.fao.org/faostat/en/#data/ET
2. Data Understanding¶
2.1 Data Understanding - Emmisions data (CEDS)¶
# load emmisions data
df_emmisions = pd.read_csv(base_path+"/air-pollution_2019.csv", sep=',')
df_emmisions.sample(5)
df_emmisions.info()
print("shape:", df_emmisions.shape)
df_emmisions.describe()
Nitrogen oxides (NOx) - gases consisting of nitrogen and oxygen, cause smog and acid rain, are emitted by vehicles and industrial processes.
Sulphur dioxide (SO2) - a colourless gas, causes acid rain, emitted from burning fossil fuels, especially coal.
Carbon monoxide (CO) - a colourless and odourless gas, toxic to humans, emitted from vehicles and the incomplete combustion of fossil fuels.
Organic carbon (OC) - a form of carbon, major component of particulate matter, emitted from vehicles and industrial processes.
NMVOCs (Non-Methane Volatile Organic Compounds) - a group of organic gases that contribute to the formation of smog and particulate matter and are emitted by industrial processes, vehicles and solvents.
Black carbon (BC) - a type of particulate matter composed mainly of carbon and produced by the incomplete combustion of fossil fuels, biofuels and biomass.
Ammonia (NH3) - colourless, strong-smelling gas that contributes to particulate matter and acid rain and is released by agricultural activities, especially livestock and fertilisers, industry and waste management.
(Hoesly et al, Historical (1750–2014) anthropogenic emissions of reactive gases and aerosols from the Community Emissions Data System (CEDS). Geosci. Model Dev. 11, 369-408, 2018a.)
sn.boxplot(data=df_emmisions, orient="h", palette="Set2")
Lets use the zscore to see the scattering of the column values in the dataframe
from scipy.stats import zscore
def abs_zscore(x):
return np.abs(zscore(x))
df_numeric = df_emmisions.select_dtypes(include=np.number)
df_zscore = df_numeric.apply(abs_zscore)
threshold = 5
outliers = df_emmisions[(df_zscore>threshold).any(axis='columns')]
none_outliers = df_emmisions[~(df_zscore>threshold).any(axis='columns')]
def plot_outliers(outliers, none_outliers, number_of_cols, number_of_rows, figsize=(24, 10)):
fig, axes = plt.subplots(number_of_cols, number_of_rows, figsize=figsize)
axe = axes.ravel()
# assign the plot to each subplot in axe
for i, c in enumerate(outliers.columns):
outliers[[c, 'Year']].plot(kind='scatter', x=c, y='Year', ax=axe[i], color='red')
none_outliers[[c, 'Year']].plot(kind='scatter', x=c, y='Year', ax=axe[i])
if c == 'Entity':
axe[i].set(xticklabels=[])
plot_outliers(outliers, none_outliers, 3, 3, (24, 10))
# get outlier values of Entity
unique_entities = outliers['Entity'].unique()
unique_entities, outliers[['Carbon monoxide (CO)', 'Year', 'Entity']].describe()
As we can see, the dataset also contains calculations for complete continents and grouped countries.
We only want to research on a country level. Thats why these values can be removed.
none_countries = ['Asia', 'Africa', 'Europe', 'High-income countries',
'Lower-middle-income countries', 'North America', 'South America',
'Upper-middle-income countries', 'World']
filtered_outliers = outliers[~outliers['Entity'].isin(none_countries)]
filtered_outliers[['Carbon monoxide (CO)', 'Year', 'Entity']].describe()
import plotly
plotly.offline.init_notebook_mode()
plotly.express.choropleth(data_frame=df_emmisions[~df_emmisions['Entity'].isin(none_countries)],
locationmode='country names',locations='Entity', color='Carbon monoxide (CO)',animation_frame='Year',
title='Carbon emission by country and year',color_continuous_scale=plotly.express.colors.sequential.Plasma,range_color=(100, 0))
fig, axes = plt.subplots(1, 2, figsize=(22, 7))
df_emmisions[df_emmisions['Entity']=='Germany'].plot(x='Year', ylabel="Emissions in million tons", title='Emissions in Germany', ax=axes[0])
df_emmisions[df_emmisions['Entity']=='United States'].plot(x='Year', ylabel="Emissions in million tons", title='Emissions in the USA', ax=axes[1])
2.2 Data Understanding - Temperature data (FAOSTAT)¶
# load emmisions data
df_temperature = pd.read_csv(base_path+"/temperature_data_2019.csv", sep=',')
df_temperature.sample(5)
df_temperature.info()
print("shape:", df_temperature.shape)
def temperature_df_preperation(df):
df = df.drop(df[df_temperature.Element == "Standard Deviation"].index)
df = df.rename(columns={'Value': 'Temperature Change'})
return df
df_temperature = temperature_df_preperation(df_temperature)
df_temperature[df_temperature['Area']=='Germany'].plot(x='Year', y="Temperature Change", ylabel="°C", figsize=(10, 5), title='Temperature Change in Germany')
df_temperature_numeric_only = df_temperature.select_dtypes(include=np.number)
df_temperature_zscore = df_temperature_numeric_only.apply(abs_zscore)
threshold = 5
outliers_temperature = df_temperature[(df_temperature_zscore>threshold).any(axis='columns')]
none_outliers_temperature = df_temperature[~(df_temperature_zscore>threshold).any(axis='columns')]
plot_outliers(outliers_temperature[["Year", "Temperature Change"]], none_outliers_temperature[["Year", "Temperature Change"]], 1, 2, (20, 6))
plotly.express.choropleth(data_frame=df_temperature, locationmode='country names', locations='Area',
color='Temperature Change', animation_frame='Year', title='Annual temperature change per country',
color_continuous_scale=plotly.express.colors.sequential.Plasma, range_color=(100, 0))
3. Data Preparation¶
We need to clean columns and rows that are not intresting in both datasets.
Additionally we want to merge both dataframes for further analysis.
# check for null values in both datasets
print(df_emmisions.isnull().sum())
print(df_temperature[["Year", "Temperature Change", "Area"]].isnull().sum())
# drop null values
df_temperature = df_temperature.dropna()
We want to merge the two dataframes by Year->Year and Entity->Area unfortunaly it looks like not all Entity labels are the same as the Area labels.
In the following cell you can see the non-existing labels in the temperature dataframe, that are missing in the emissions dataframe. A "difflib" function was used to see the closest matches of missing labels. With this result the temperature dataframe was manually edited to use the same labels as the emission dataframe. But as you can see, not all labels could be found.
# compare Area and Entity in both data sets
import difflib
df_emmisions_countries = df_emmisions.Entity.unique()
df_temperature_countries = df_temperature.Area.unique()
diff_countries = set(df_temperature_countries).difference(set(df_emmisions_countries))
for country in diff_countries:
print(country, difflib.get_close_matches(country, df_emmisions_countries))
# read in edited temperature dataframe
df_temperature = pd.read_csv(base_path+"/temperature_data_2019_country_corrected.csv", sep=',')
df_temperature = temperature_df_preperation(df_temperature)
df_temperature = df_temperature.dropna()
# merge data sets
new_df = pd.merge(df_emmisions, df_temperature, how='inner', left_on=['Entity','Year'], right_on = ['Area','Year'])
# remove unnecessary columns
new_merged_df = new_df[['Nitrogen oxide (NOx)', 'Sulphur dioxide (SO₂)', 'Carbon monoxide (CO)','Organic carbon (OC)', 'NMVOCs', 'Black carbon (BC)', 'Ammonia (NH₃)', 'Year', 'Temperature Change', 'Entity']]
new_df.columns, new_merged_df.columns
Let's create a simple linear correlation heatmap. Between the pollution values and the temperature change value. The data is currently not grouped by country.
plt.figure(figsize=(8, 6))
new_df_cols = new_df[['Nitrogen oxide (NOx)', 'Sulphur dioxide (SO₂)', 'Carbon monoxide (CO)','Organic carbon (OC)', 'NMVOCs', 'Black carbon (BC)', 'Ammonia (NH₃)', 'Year', 'Temperature Change', 'Entity']]
heatmap = sn.heatmap(new_df_cols.corr()[['Temperature Change']].sort_values(by='Temperature Change', ascending=False), vmin=-1, vmax=1, annot=True, cmap='BrBG')
heatmap.set_title('Features Correlating with temperature change', fontdict={'fontsize':18}, pad=16);
Let's try this correlation analysis with data that comes only from the country the emission were produced in.
# create a correlation heatmap for the individual countries
corrs_per_country = new_df_cols.groupby(['Entity']).corrwith(new_df_cols['Temperature Change'])
heatmap = sn.heatmap(corrs_per_country.iloc[65:73].sort_values(by='Temperature Change', ascending=False), vmin=-1, vmax=1, annot=True, cmap='BrBG')
# define final dataframe for data modelling
df_prepared = new_df[['Nitrogen oxide (NOx)', 'Sulphur dioxide (SO₂)', 'Carbon monoxide (CO)','Organic carbon (OC)', 'NMVOCs', 'Black carbon (BC)', 'Ammonia (NH₃)', 'Year', 'Entity', 'Temperature Change']]
# prepare a model (58 Years of data per country) and leave out 2019 for later evaluation
df_regression_germany = df_prepared[(df_prepared['Entity']=='Germany') & (df_prepared['Year']<2019)].sort_values(['Entity', 'Year'])
df_regression_all = df_prepared[(df_prepared['Year']<2019)].sort_values(['Entity', 'Year'])
df_evaluation_germany = df_prepared[(df_prepared['Entity']=='Germany') & (df_prepared['Year']>=2019)].sort_values(['Entity', 'Year'])
df_evaluation_all = df_prepared[(df_prepared['Year']>=2019)].sort_values(['Entity', 'Year'])
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
scaler = MinMaxScaler()
df_regression_x = df_regression_germany[['Nitrogen oxide (NOx)', 'Sulphur dioxide (SO₂)', 'Carbon monoxide (CO)','Organic carbon (OC)', 'NMVOCs', 'Black carbon (BC)', 'Ammonia (NH₃)', 'Year']]
df_regression_y = df_regression_germany[['Temperature Change']]
X_train, X_test, y_train, y_test = train_test_split(df_regression_x, df_regression_y, test_size=0.2, random_state=42, shuffle=True)
# later used for evaluation
X_test_original = X_test.copy()
# min max scale data between 0 and 1
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
y_train = scaler.fit_transform(y_train)
y_test = scaler.transform(y_test)
4. Modeling¶
Use a machine learning "Random Forest" model for regression.
Tune Hyperparameters with random cv search. (reaches gradient descent for best values faster)
Deep Learning - Bengio, Goodfellow, Courville
For the regression a simple RandomForestRegressor from the sklearn lib was used.
https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html
A Random Forest is a type of machine learning algorithm that builds multiple decision trees (hence "forest") and combines their predictions to classify new data points. Each decision tree is built using a random sample of the data and a random subset of the features, which helps to reduce overfitting and improve the overall accuracy of the model. The final prediction is made by taking a majority vote of the predictions from all the decision trees in the forest. It is called random because it randomly selects observations and features to build decision trees.
https://towardsdatascience.com/understanding-random-forest-58381e0602d2
# random search for RandomForestRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import RandomizedSearchCV
from pprint import pprint
# Number of trees in random forest
n_estimators = [int(x) for x in np.linspace(start = 4, stop = 50, num = 4)]
# Number of features to consider at every split
max_features = ['auto', 'sqrt']
# Maximum number of levels in tree
max_depth = [int(x) for x in np.linspace(2, 10, num = 2)]
max_depth.append(None)
# Minimum number of samples required to split a node
min_samples_split = [2, 4, 8]
# Minimum number of samples required at each leaf node
min_samples_leaf = [1, 2, 4]
# Method of selecting samples for training each tree
bootstrap = [True, False]
# Create the random grid
random_grid = {'n_estimators': n_estimators,
'max_features': max_features,
'max_depth': max_depth,
'min_samples_split': min_samples_split,
'min_samples_leaf': min_samples_leaf,
'bootstrap': bootstrap}
# First create the base model to tune
rf = RandomForestRegressor()
# Use the random grid to search for best hyperparameters
rf_random = RandomizedSearchCV(estimator = rf, param_distributions = random_grid, n_iter = 3, cv = 3, verbose=0, random_state=42, n_jobs = -1)
# Fit the random search model (ravel is used to flatten the y_train array)
rf_random.fit(X_train, y_train.ravel())
best_params = rf_random.best_params_
print(best_params)
# RandomForrestRegressor
print(best_params)
#regressor = RandomForestRegressor(n_estimators=best_params['n_estimators'], min_samples_split=best_params['min_samples_split'], min_samples_leaf=best_params['min_samples_leaf'], max_features=best_params['max_features'], max_depth=best_params['max_depth'], bootstrap=best_params['bootstrap'])
# actual model that will be using for the evaluation
regressor = RandomForestRegressor(n_estimators=4, min_samples_split=2, min_samples_leaf=4, max_features="sqrt", max_depth=10, bootstrap=False)
regressor.fit(X_train, y_train.ravel())
feature_importances = pd.DataFrame(regressor.feature_importances_,
index = df_regression_x.columns,
columns=['importance']).sort_values('importance', ascending=False)
feature_importances.plot(kind='barh', figsize=(10, 5), title='Feature importance')
5. Evaluation¶
# predict values
y_pred = regressor.predict(X_test)
# rescale results and reshape to 2d array
y_pred = scaler.inverse_transform(y_pred.reshape(-1, 1))
y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))
X_test_rescaled = scaler.inverse_transform(X_test.reshape(-1, 1))
print('Mean Squared Error:', mean_squared_error(y_test, y_pred))
# plot results
plt.figure(figsize=(10, 5))
plt.plot(y_test_rescaled, color = 'red', label = 'Real Temperature Change')
plt.plot(y_pred, color = 'blue', label = 'Predicted Temperature Change')
plt.title('Temperature Change Prediction')
plt.xlabel('Sample Count')
plt.ylabel('Temperature Change')
plt.legend()
plt.show()
# add Temperature Change Prediction to X_test_original
X_test_original['Temperature Change Prediction'] = y_pred
X_test_original['Temperature Change'] = y_test_rescaled
print(X_test_original[['Year', 'Temperature Change', 'Temperature Change Prediction']])
def model_by_country(df):
df_regression_x = df[['Nitrogen oxide (NOx)', 'Sulphur dioxide (SO₂)', 'Carbon monoxide (CO)','Organic carbon (OC)', 'NMVOCs', 'Black carbon (BC)', 'Ammonia (NH₃)', 'Year']]
df_regression_y = df[['Temperature Change']]
X_train, X_test, y_train, y_test = train_test_split(df_regression_x, df_regression_y, test_size=0.2, random_state=42, shuffle=True)
regressor = RandomForestRegressor(n_estimators=4, min_samples_split=2, min_samples_leaf=4, max_features="sqrt", max_depth=10, bootstrap=False)
return regressor.fit(X_train, y_train.values.ravel())
model = df_prepared.sort_values(['Entity', 'Year']).groupby(['Entity']).apply(model_by_country)
model = model.reset_index()
model = model.rename(columns={0: 'Temperature Change Model'})
results_2019 = []
for index, row in df_evaluation_all.iterrows():
row_x = [row['Nitrogen oxide (NOx)'], row['Sulphur dioxide (SO₂)'], row['Carbon monoxide (CO)'], row['Organic carbon (OC)'], row['NMVOCs'], row['Black carbon (BC)'], row['Ammonia (NH₃)'], row['Year']]
row_y = row['Temperature Change']
model_country = model[(model['Entity']==row['Entity'])]
prediction = model_country.iloc[0]["Temperature Change Model"].predict([row_x])
difference = np.abs(prediction - row_y)
error = difference / prediction
mean_squared_error = np.mean(error)
results_2019.append([row['Entity'], row['Year'], prediction[0], row_y, mean_squared_error])
# results_2019 to dataframe
df_results_2019 = pd.DataFrame(results_2019, columns=['Entity', 'Year', 'Prediction', 'Real', 'Mean Squared Error'])
# df_results_2019 avg for Mean Squared Error
print("Error across all countries:",df_results_2019['Mean Squared Error'].mean())
# plot results of 2019 seaborn
plt.figure(figsize=(10, 5))
sn.set_style("whitegrid")
sn.barplot(x="Entity", y="Mean Squared Error", data=df_results_2019[68:78])
plt.xticks(rotation=90)
plt.title('Mean Squared Error of 2019')
plt.xlabel('Country')
plt.ylabel('Mean Squared Error')
plt.show()
So what does this result mean for our original question "We see that the original question "Does the air pollution of a country have a direct influence on the annual average temperature?" ?
We can see that it is possible for some countries better for others worse to use the emission values of the current year to predict the temperature change of the country. That means emissions or air pollution can or have a effect on our climate/temperature.
6. Deployment¶
Now we would like to make our country based prediction models available for scientist. So that they can use measured emission values of their country to predict the temperature change for a year.
Therefore we will use flask a simple python web framework to make this possible.
https://flask.palletsprojects.com/en/2.2.x/